clear all

global Dividend

Dividend=0;

% Given Par2 optimal:
n = 50; tau = 0.1;
r = 0.036;
Par1 = zeros(25,1);
load ParEst_eh0_time0_Mex ParEst_eh0_time0_Mex
load ParEst_eh0_time1_Mex ParEst_eh0_time1_Mex
Par1(1:24) = ParEst_eh0_time0_Mex;
Par1(21) = 10^Par1(21);
Par1(22:24) = Par1(22:24)*1000;
Par1(25) = 0;
% simulation here
Par1(20)=Par1(20)*1.1;
gamma =  Par1(25); a = Par1(22); b1 = Par1(23); b2 = Par1(24);
alphaf = Par1(17); betaf = Par1(18); alphai = Par1(19); betai = Par1(20); theta = Par1(21);
load W_time1_age0_eh0_comp11.out
w_f_g=W_time1_age0_eh0_comp11; % Read wages, g and f distributions
wf_1=w_f_g(:,1);        % Wage in formal sector
wi_1=w_f_g(:,2);        % Wage in informal sector
wf_2=w_f_g(:,7);      % Wage in formal sector
wi_2=w_f_g(:,8);      % Wage in informal sector

[Ff_1,ff_1,Fi_1,fi_1,Ff_2,ff_2,Fi_2,fi_2]=solve_wage_offer_dist(Par1,wf_1,wi_1,wf_2,wi_2,n);
[Wff,Wfi,Wfn,Wif,Wnf,Wii,Win,Wni,Wnn, Rw_fi_ni, Rw_ii_ni, Rw_ii_in, Rw_if_in,...
    Rw_fi_fn, Rw_ff_fn, Rw_if_nf, Rw_ff_nf, Rw_in_nn, Rw_fn_nn, Rw_ni_nn, Rw_nf_nn...
    ]=solve_value_function_A3_C(Par1,Ff_1,wf_1,Fi_1,wi_1,Ff_2,wf_2,Fi_2,wi_2,n);
[HHstates,wage_moments,transitions,gff, gfi, gif, gii, gfn, gnf, gin, gni]=data_simulation_A12_end(Par1,ff_1,wf_1,fi_1,wi_1,ff_2,wf_2,fi_2,wi_2,Wff,Wfi,Wfn,Wif,Wnf,Wii,Win,Wni,Wnn);
wage_moments_beta_i_plus10 = wage_moments;
wage_moments_beta_i_plus10_red=wage_moments(1:8);%meanwi_1; meanwf_1; meanwi_2 ; meanwf_2; sdwi_1; sdwf_1; sdwi_2; sdwf_2
HH_states_beta_i_plus10 = HHstates;
save wage_moments_beta_i_plus10_red wage_moments_beta_i_plus10_red
save HH_states_beta_i_plus10 HH_states_beta_i_plus10

mff=HHstates(1);
mfi=HHstates(2);
mfn=HHstates(3);
mif=HHstates(4);
mnf=HHstates(5);
mii=HHstates(6);
min_=HHstates(7);
mni=HHstates(8);
mnn=HHstates(9);
Un_Head = mnf + mni + mnn; Un_Spouse = mnn + min_ + mfn; Inf_rate = (mnn + min_ + mni + mii);
Agg_stocks_beta_i_plus10= [Inf_rate; Un_Head; Un_Spouse];
save Agg_stocks_beta_i_plus10 Agg_stocks_beta_i_plus10

% Transitions 
Trans_head = [transitions(3); transitions(4); transitions(1); transitions(5); transitions(2); transitions(6); transitions(13); transitions(14)];
Trans_spouse = [transitions(9); transitions(10); transitions(7); transitions(11); transitions(8); transitions(12); transitions(15); transitions(16)];
Transitions_beta_i_plus10 = [Trans_head;NaN;NaN;Trans_spouse];
save Transitions_beta_i_plus10 Transitions_beta_i_plus10


